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The Onsager reaction-field (ORF) theory is extended so 
that it applies to any three-dimensional Bravais lattice with 
a basis. The ORF calculation is used to predict the critical 
temperature for classical Ising, XY, and Heisenberg magnetic 
models, in particular, on diamond and hexagonal close packed 
lattices. Results are compared with series extrapolations and 
other theoretical approaches where available. For the hep 
lattice the ORF calculation is seen to be exactly equivalent 
to a Green's function approach. 

PACS numbers: 75.10.Hk, 75.10.-b, 75.40. Cx 



I. INTRODUCTION 

The Onsager reaction field (ORF) theoryB is an im- 
proved form of mean-field theory that includes at least 
partially the effects of correlations between nearby atoms, 
ft was Qciginally applied in magnetism by Bcout and 
Thomas,0 and more-, recently to spin r, glasses itiner- 
ant electron. systemsB Hubbard models,0 and anisotropic 
fleisenbergQ and XYO models. The procedure is versatile 
and has been used to estimate specific heat, susceptibil- 
ity and correlations above the critical temperature, Tc, 
as well as Tc itself. A clear review of the ORF method 
applied ta|-|three-dimensional (3D) Ising models is given 
by WhiteJa The method has been applied on the stan- 
dard Bravais lattices, including simple cubic (sc), body- 
centered-cubic (bcc), and face-centered-cubic (fee), with 
results given in terms of integrals over the associated Bril- 
louin zone (BZ). However, a modification of these calcu- 
lations is needed to consider other lattices which are not 
in the Bravais classification. Here we show how to ap- 
ply the ORF procedure to any non-Bravais lattice that 
can be considered as an underlying Bravais lattice with 
a basis. In particular, the diamond and hexagonal close 
packed lattices are analyzed, both of which have two- 
atom bases. I-I 

In the usual mean-field theory due to Weiss,cl a chosen 
atom (or spin, for the magnetic problems we consider) 
is viewed as interacting with the average, or mean-field, 
of its nearest neighbors. The exact Hamiltonian is re- 
placed by the mean-field one, in which the neighbors are 
introduced as the mean-field acting on the central atom. 
However, the central atom itself influences the neighbors, 
and therefore the mean-field usually includes a part di- 
rectly attributed to the central atom. This means the 
mean-field includes a part that might be considered a 
self-interaction effect, which should really be subtracted 



out. This results in an overestimate of the critical tem- 
perature, Tc- The ORF procedure is simply a way to 
estimate and subtract out this self-interaction part, i.e., 
by adding a "reaction field" term that accomplishes this. 
In this way, the estimate of Tc is brought down, indeed, 
usually ORF leads to an underestimate of Tc- 

The standard ORF approach uses as input the specific 
lattice structure, be it sc, fee, bcc, etc. However, the 
usual approach and well-known formulas require the per- 
fect periodicity of a Bravais lattice-all atoms are taken 
as equivalent. On the other hand, we haMe been inter- 
ested in mean-field and other calculationsEJ for diamond 
lattices because of the low coordination number (z = 4). 
The diamond lattice is not a Bravais lattice-all sites do 
not have the same surroundings, instead, the diamond 
lattice can be considered to be a fee lattice with a two 
atom basis. Thus it is interesting to understand how to 
apply the ORF procedure to such a system. Recently 
there is interest in fecromagnetic ordering of hep '^He 
at low temperatures,Eil assumed to be described by a 
Heisenberg model. The hep lattice is another example 
of a non-Bravais lattice; it can be considered as simple 
hexagonal (stacked triangular nets) also with a two-atom 
basis. Here we show how to extend the standard ORF 
calculation of Tc to these two systems, however, our ap- 
proach will apply to any Bravais lattice with a basis, i.e., 
any system with multiple atoms per unit cell. 

For simplicity we display formulas for Ising models on 
a 3D lattice with spins Sn — ±5" and coordination num- 
ber z. However, the modifications to treat n-component 
spins [i.e., XY (n=2), Heisenberg (n=3), etc. ] are min- 
imal and will be noted where they are appropriate. The 
Hamiltonian is 



(1) 



where each sum is over all of the lattice sites, and the 
bond coupling strength Jn.m depends only on the neigh- 
bor displacement, n — m, and is of the same strength 
J for all near neighbor pairs. iJ„ is a spatially varying 
applied field. 



II. ONSAGER REACTION FIELD CORRECTION: 
BRAVAIS LATTICES 

In the ORF calculation (See Ref. ^for more details), a 
spin at a chosen site interacts with the mean-field reduced 
by "reaction field" that depends on the spin at that 
site.EJ For completeness we summarize key aspects of this 
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calculation here to see how the extension to a non-Bravais 
lattice is accomplished. 

To effect the reaction term, in the real space Hamilto- 
nian an extra self-interaction term is added: 



A ^ ^ SyiSn 



(2) 



This is equivalently a delta-function exchange term of 
strength A. The constant A is the reaction-field, which 
is determined self-consistently in the calculation, by a 
constraint on the magnetic susceptibility, below. 

The Hamiltonian Ti. = Tio + Tirf, transformed into 
wavevector space is. 
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H = - - [5_q(Jq - A)5q + + (3) 

q 

where the Fourier-space quantities derive from 
1 



5„ = 



Jn,m — Jq 



,iq-(n-m) 



(4) 



(5) 



together with a similar definition for the field = i?_q, 
for real applied field Hn- N is the number of lattice sites. 
The sum in Eq. ^ is over all q in the appropriate Brilluoin 
zone. The Fourier-transformed exchange interaction is 



Jq — Jr 



(6) 



i.e., a sum over displacements to nearest neighbors, r = 
m — n. 

The q-dependent magnetization and zero-field suscep- 
tibility have the usual definitions. 



Mq = (5q 



9Mq 



PiS-qSq) 



(7a) 



(7b) 



Then for an arbitrary Hamiltonian, there results the con- 
straint (used to determine the reaction field A) 



n 



Xo 



(8) 



The number of spin components n enters here, in the ex- 
pression for xoi when models other than the Ising model 
are considered. 

Now we consider the mean-field approximation in 
q-space, i.e., using the random phase approximation 
(RPA), magnetization components at different wavevec- 
tors are assumed to be independent. The mean- 
field Hamiltonian for the interaction of a negative q- 
component can be written 



n 



(9) 



where the effective (or mean-field) magnetic field is 

Hf = H^ + (Jq - A) (Sq) (10) 



From the definition (7b), and the relation. 



the mean-field Hamiltonian gives the well-known expres- 
sion, 



Xq 



Xo 



1 - Xo(^q - A) 



(12) 



The reaction field A is determined by forcing this expres- 
sion for Xq to satisfy the constraint (^) . The critical tem- 
perature Tc is determined as the temperature at which 
Xq=o diverges. For lower temperatures the ORF calcu- 
lation gives a negative susceptibility at q = 0, signifying 
the presence of the ordered state. In a continuum limit 
of the constraint (^), a short manipulation leads to an 
expression for the critical temperature. 



JoS^ 
nl 



(13) 



The constant / is given from a q-space integral over the 
appropriate Brilluoin zone: 



Jo 



BZ 



(2^)3 Jo-Jq 



(14) 



and V/N is the specific volume per lattice site (for ex- 
ample, V/N = a^, ifl^, ja^ for sc, bcc and fee lattices, 
respectively, where is the cubic unit cell volume). 
At q = we have Jq = zJ, which gives the energy 
scale in the mean-field approximation. Then the di- 
mensionless integral / gives the correction to the mean- 
field critical temperature. From integration over the 
appropriate Brilluoin zones, the values of / are 1.517, 
1.393 and 1.343 for sc, bcc and fee lattices, respectively. 
When applied to the 3D Ising model (ti = 1) one getsB 
keTc/JS^ = z/I = 3.955,5.743,8.932 for sc, bcc and 
fee lattices, considerable improvements over the standard 
mean-field results given from ksTc/JS^ = z. They com- , 
pare favorably with the exact Ising model resultal3^t2l 
from series: ksT^/JS'^ = 4.5103,6.3508,9.794, respec- 
tively. For the Heisenberg model [n = 3), the ORF 
predictions would be ksTc/JS"^ = 1.318, 1.914, for sc 
and bcc lattices, whereas precise Monte Carlo estimatesEil 
give ksTc/JS^ = 1.443,2.054, respectively. Note, how- 
ever, that the ORF estimates are all below the exact 
results. 
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III. BRAVAIS LATTICES WITH A BASIS 

Now suppose the lattice has an underlying set of N 
Bravais lattice points n, each of which has a two-atom 
basis {0, d}. (The generalization to a larger basis is 
straightforward.) Therefore at each site n we suppose 
there is a two component field written as a column 
vector: 



Sn+d 



(15) 



It is useful to employ a sublattice notation, Sn = 
Sn,Sn+d = Sn ■ The exchange interaction occurs be- 
tween neighboring Bravais sites, via a 2 x 2 matrix, Gn,m- 



Ti-cx — - ^ ^ ^ ■ Gn,m ■ 



(16) 



where it is stressed that the sums are over all Bravais 
sites; the factor of 1/2 avoids double counting, and 



qBA qBB 



<Ai,m <-^n,m+d 
</n+d,m <Ai+d,m+d 



(17) 



In fact, the matrix Gn.m is taken as zero unless n — m 
is a near neighbor displacement. The details of the spe- 
cific lattice will determine which components of G are 
nonzero. 

The reaction terms, two for each site n, can be written 
with a 2 X 2 unit matrix. 



"Worf — 



1 
1 



(18) 



which is equivalent to expression (g), and essentially 
shifts the original exchange matrix by — 2A(5„.mX, where 
I is the 2x2 unit matrix. 

Finally, for the purpose of the calculation, there is a 
separate applied field for each sublattice, so at a given 
site n, we have fields (acting on S'„) and (acting 
on S'n+d)- These compose a column vector field, 



(19) 



Then we choose to write the q-space Hamiltonian, in- 
cluding the reaction term and applied fields, as 



q 



(20) 



The Fourier transforms needed here obviously are related 
to those already defined for J„_in, 'S'n7 etc. The most sig- 
nificant difference from Eq. (||) is the presence of "2" on 



A, due to the two-atom basis. This Hamiltonian is ex- 
act. Now we define the sublattice magnetizations (where 

Ml, = {S\) (21) 
and related susceptibilities (where also j = A,B), 



(22) 



In the RPA, from the point of view of the negative 
q-components, the Hamiltonian is approximated as 

^ = - E [^"^q • ^ + ^-q • ^q • <^q)] (23) 



where we use the shorthand notation for the shifted ex- 
change interaction. 



Kq = Gq- 2X1 



(24) 



The above Hamiltonian can alternatively be written in 
terms of effective fields in components. 



(25) 



< = < + A--(5-) + <«(5«> (26a) 



^AA I oA 



AB I oB\ 



H. 



B cff 



(26b) 



Then in the limit of zero applied field, using this RPA 
Hamiltonian, the susceptibility definitions (E2[) become 



dM' dH„ 

— q_ SL 

Aq 



icff 



(27) 



and we get equations for the susceptibility components, 

X^^ = P{S^S\), [1 + K^^x^^ + K^^x^^^] (28a) 



(28b) 



xr - P{S^S^^)o [1 + Kl^xT + ^''xr ] (28c) 



^ = f3{S^S^^)o [<^xf ^ + K^^i^] (28d) 

where ( )o means the expectation value using the RPA 
Hamiltonian. In the high-temperature limit, these ex- 
pectations are 

P{S^S\)o - PiS^S^^o = -S' ^ xo- (29) 
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The equations ( pq ) can be solved in the general case for 
all four susceptibility components. We get 



Xq^ = XO 



1 



BA 
Aq 



Xq 



j^BA AA 
1 - XOKBB 



(30a) 



(30b) 



and similar equations for Xq ^-nd Xq by appropriately 
interchanging the AB indeces. For the lattices considered 
in this paper, however, there are the symmetries, = 



K^^ and K^^ 



Therefore the solutions are 



seen to satisfy symmetries = Xa^ i Xq^ — (Xq^)* 



Now consider how to determine the reaction field A, 
and subsequently, Tc- A little consideration using the 
definitions of Xq^ ^^id shows that there is still the 
constraint. 



-Yx^^ 



Xo 



(31) 



Converting to the continuum limit, 
3 ^ 

' d-'q Xq^ = Xo 



(32) 



Jo — Gq 



r'ABr'BA 



(35) 



The integral / defined in this way again gives the cor- 
rection to the mean- field prediction for Tc. Thus the 
determination of Tc has been reduced to evaluating this 
integral over the Brilluoin zone. 

Below we will use Eq. (|4|) to estimate the critical tem- 
perature for diamond and hep lattices. But first we can 
verify that the result is correct by using it for a bcc lat- 
tice, considered as a simple-cubic lattice with a basis, 
where we already know the standard ORF result for Tc. 



IV. BCC LATTICE AS SC WITH BASIS 

The bcc lattice points can be generated from the simple 
cubic primitive vectors, ai = ax, sl2 — ay, and 33 = 
az, where a is the cubic cell lattice constant, together 
with the basis, {0,d}, where d = |-(i + y + z) is the 
displacement to the body-centered point. The lattice can 
be thought of as a pair of interpenetrating sc lattices (A, 
B sublattices) with displacement d. The near neighbors 
of an 'A' site are all "B' sites, and vice-versa, with the 
result that the G^^ and G^^ couplings are all zero. The 
nonzero ^ = G^^ couplings depend only on the near 
neighbor displacement, r = m — n, as follows: 

Gf = J, r = 0, ai, a2, as, ai + a2 + as, 

ai -I- a2, a2 4- as, as + ai (36a) 



This equation implicitly determines the reaction coupling 
A for any T > Tc. It is not clear how to get an explicit 
solution for A from it; a solution for A(T) for a given 
lattice can be found numerically (below). 

Now, just as described for the Bravais lattice systems, 
the critical temperature Tc is the temperature at which 
any of the susceptibility components, at q = 0, diverges. 
(The Xq are well defined on the high temperature side 
of Tc.) So Tc is determined as the point at which the de- 
nominator of Eq. (30a) goes to zero, leading to a relation 
between the critical temperature (via xo) and the critical 
reaction field, 



Xo ^ + 2A - G^^ 



r'ABr'BA 



(33) 



Using this in th e constraint equation ( P2|) together with 



the result (|30a|) for Xa^^ gives the general result when 



the symmetry G^^ = G^^ holds. 



G 



AB 



J, r 



0, -ai, 
ai - a2. 



a2, -as, 
- a2 - as. 



- ai - a2 - as, 
as - ai (36b) 



where the terms for r = correspond to the coupling 
within the two-atom basis. The fact that G^ m — 
G^^ = simplifies the determination of Tc consider- 
ably, as we only need to know the product, G^^G^^. 

The Fourier-transformed interactions are found from 
sums over all the nonzero Gr (15 possible terms): 



G, 



AB 



ABJqr 



(37) 



G 



AB 



J{1 



^g-*(<?^+9i/) _|_ g-^iqy+q^) ^_ g~Hq^+qm) 
+e-^(9.+'3«+9.)} (38) 



knTc = 



nl 



Ja{Jo — G„ ) 



N 7bz (2^)3 (Jo - G^^r ~ G^BQBA 
where Jq is the effective q = exchange strength. 



(34a) 



(34b) 



In this and the following equations, qx,<ly, 9z are in units 
of 1 /a. For the underlying sc lattice, the density of points 
is 1 for every volume a'^, i.e., V/N — a^. After a short 
calculation, there results 



q — 8J {1 -f cos(72, -f cos(7y -I- cosg^ 



qABqBA 

cos qx cos 
- cos qx cos 



+ cos qy cos qz 4- cos qz cos qx 
cosg^} (39) 
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and the determination of Tc 
simplified integral, 



relies on evaluation of the 



/ = 



rf3 



1 



■^0 



(27r)3 J2_G^SGBA 



(40) 



where Jn 



— ^3 



'-'0 



8J = zJ and the factor 



was absorbed into the dimensionless q. 

" dqx_ r dqy_ r dqz_ 

1 1 — o [l + cos qx + cos qy + cos qz 
L 8 

+ cos qx cos qy + cos qy cos qz + cos qz cos qx 



- cos qx cos gy cos qz] | 



(41) 



The integral / was evaluated numerically by sampling 
qx,qy,qz uniformly with a constant increment, and then 
using an extrapolation of the results in the limit that the 
increment goes to zero. Essentially, we let dqx — dqy — 
dqz = 2Tr/Nx, where Nx is some integer number of di- 
visions of the axes, and then generated a simple cubic 
lattice of sampling points from these. The integral / was 
then estimated as a sum over the resulting cubic grid of 
sampling points in the BZ in q-space. A plot of / ver- 
sus \/Nx results in a straight line whose extrapolation 
to \/Nx gives a very accurate estimate of the inte- 
gral. (Errors in the integral estimate clearly go as 1/Nx-) 
In this way we found / = 1.39321. Therefore, the esti- 
mate of critical temperature for the bcc system obtained 
viewing it as sc with a basis is 



8JS^ 



1 



5.74213 



1.39321 



(42) 



For comparison with the standard ORF procedure, 
when the original bcc lattice is used, the specific vol- 
ume is V/N — ifl^. The Fourier transformed exchange 
(Eq. I) is 



r o T <lx Qy qz 
Ja = oJ cos COS — COS 

^ 2 2 2 



(43) 



The integral needed for Eq. (p4) is 



f ^'g r 

Jbz {27r)A 



qx 



qy 



r 



= - X 2.786 = 1.393, 



(44) 



where the integral is over the Brilluoin zone for the bcc 
lattice (an fee Wigner-Seitz cell), and was evaluated by 
the sampling technique already described. The integral is 
exactly twice that given by expression (^), but because 
V/N = is half as large for the bcc lattice as for the 
sc lattice, the result for from Eq. (^) is that given by 
Eq. (^). Thus this approach using the two-atom basis, 
for the bcc system, is equivalent to the standard ORF 



procedure, and should be reliable for application to the 
diamond and hep lattices. 

As an interesting mathematical note, the integral / in 
(^) can be rewritten using the periodicity and symmetry 
in q-space to be over the same region as in ( ^ ) . A cubic 
ceh -2tt < qx < 2tt, -27r < qz < 27r, -27r < qz < 27r, 
contains 4 copies of the Brilluoin zone of the bcc lattice. 
Thus we can integrate over this cubic cell and divide by 
4; also, shifting the variables of integration to q'x — qx/2, 
etc., leads to an additional factor of 2^ and gives 



. 1 1 
/ = - X - X 2^ 
2 4 



^ dq;x 
27r 



2ti 



27r 



1 1 — cos q'x cos q'y cos q'^ | 



(45) 



Overall, the prefactor on the integral is one. Thus the 
integral here must be the same as the integral in (|4l]). 
The integrands, however, are not equivalent; there is no 
way to transform one into the other. 

Finally we also note that it is easier and more precise to 
evaluate ( ^5| ) by a uniform cubic grid of sampling points 
than expression (^), because a cubic grid of sampling 
points easily fits to the cubic integration boundaries. On 
the other hand, for expression ( ^ ) , there is always more 
difficulty to fit any grid of points smoothly to all of the 
integration boundaries of the fee Wigner-Seitz cell in q- 
space, leading to greater discretization errors. Removing 
these errors is important for producing a smooth extrap- 
olation to Nx — > c». 



V. DIAMOND LATTICE 

The diamond lattice can be considered as an fee lattice 
with a two atom basis. The fee primitive vectors are 
ai = |(x -I- y),a2 = §(?; + z),a2 = f (i 4- x), where a is 
the standard cubic cell lattice constant, and the basis is 
{0, d}, where d = ^{x + y + z). The nearest neighbors of 
a site n -I- d (a 'B'-site) are {n, n + ai, n-|-a2, n-f- as} (all 
'A'-sites). The nearest neighbors of a site n (an 'A'-site) 
are at {n + d, n + d — ai, n -I- d — a2, n + d — as} (all 
'B'-sites). As a result, only G^^ and G^^ are nonzero. 
In terms of the neighbor displacements r = m — n, the 
only nonzero exchange couplings are 



G^'' = J, r = {0,ai,a2,a3} 

G^^^ = J, r = {0,-ai,-a2,-a3} 
Then it is straightforward to evaluate 



G. 



r 

= j[i 



AS iq r 



{G^^r 



(46a) 
(46b) 

(47a) 
(47b) 



and what is needed for the Tc integral: 



5 



4J^ 



qx qy 

+ COS COS — 

2 2 

qy qz 

+ cos — COS — 

2 2 



cos — COS — 

2 2 



(48) 



where g's are in units of 1/a. Clearly we also have once 
again, Jq = Gq^ = Gq^ = 4J = zJ. The specific 
volume per lattice point on the underlying fee lattice is 
V/N = ia^. Then wiU be evaluated using Eq. @ 
and the simplified integral like which becomes 



4 



BZ 

f 

1-4 



(27r)3 
1 - 



qx qy 

cos — cos -7- 
2 2 



gy ?z 

+ cos — cos — 
2 2 



cos — cos 

2 2 



r 



(49) 



The BZ for the fee lattice is a bcc Wigner-Seitz cell, with 
lattice constant An /a. The easiest and most precise way 
to evaluate this integral is to apply the same procedure 
we noted for the bcc system: Change the integration re- 
gion to the cube, —2tt < Qx < 2tt, —2tt < Qy < 2tt, —2tt < 
qz < 2tt, which contains 2 copies of the BZ, therefore in- 
clude a factor of i, and sum over points on a cubic grid. 
Also make the variable change, q' = q/2, which leads to 
a factor of 2^^, giving 



4 2 ./ ^ 27r 



2-K 



2tt 



1 r 



1 + cos q'^ cos q' 



cos q' cos q'^ 



cos q'^ cos q'^ 



r 



(50) 



Using a cubic grid spacing dq'^ ~ dq'y ~ dq'^ — 2tt/Nx, 
we evaluated the integral as a sum over points within the 
cubical integration region, for ranging from 100 to 
2000. A plot of / versus gives a straight line (Fig. 

|l|), and its extrapolation to 1/N.j; gives / = 1.79288. 
Thus the critical temperature is estimated as 



1 



2.23105 



1.79288 



(51) 



The exact result for Tc (Ising model, n = 1) as esti- 
mated from series expansions, is known to be ksTc — 
2.7040JS'^. Thus the ORF calculation, as is usual, un- 
derestimates Tc but is a considerable improvement over 
the simple mean- field result, kBTc = AJS^. 



a{\x + ^y),aL'i = cz, where a and c are the lattice con- 
stants. For the hep system, a two-atom basis of {0, d} is 
used, where d = i(ai +32) + ^aa, and one triangular net 
is stacked on top of the previous one, but shifted to be 
over the centers of one set of the triangular cells below, 
in what is usually referred to as the ARAB... packing. 
For the lattice constant ratio c/a — a/8/3, the highest 
density packing is obtained, however, for the calculation 
here this ratio does not directly enter, and need not be 
specified. Instead, it is interesting to consider that the 
near neighbor exchange interactions within the triangu- 
lar nets (xy-plane) have one strength, Jxy^ while there 
is a different strength, J^, for the bonds between the 
planes. In general we can consider the calculation of Tc 
as a function of the ratio, A = J^l Jxy 

We present first the calculation of T'c(A) for the sim- 
ple hexagonal Bravais lattice, using the standard ORF 
theory in Sec. which acts as an introduction to the 
corresponding calculation for the hep system, because 
they both rely on the same information concerning the 
Brilluoin zone. 



A. Simple Hexagonal Bravais Lattice 

Here there are 6 neighbor displacements from 
some arbitrary site to neighbors in the same plane 
{±ai, ±a2, ±(ai — a2)}, with exchange strength, J^y 
The remaining two neighbors, with displacements, ±a3, 
have exchange strengths, J^. A short calculation shows 
that the q-space exchange (Eq. ||) is 

Jq ~ 2Jxy [cos q • ai + cos q • a2 + cos q • (ai — a2)] 

+2 J2 cos q • as (52a) 

1 \/3 

= 2Jxy [cos q^a + 2 cos -q^a cos — Qyo] 

-t-2J2C0sq^c (52b) 

and Jo = QJxy + 2J2 ~ 2(3 + A) J^^y will determine the 
mean-field critical temperature. 



The area of one triangle in the net is ia x 



-^a^, and there is ^-site per triangle per layer. Thus the 

specific volume per site is V/N = ^c?c. The primitive 
vectors of the reciprocal space are 



47r , V3 . 1 



(53a) 



VI. SIMPLE HEXAGONAL BRAVAIS LATTICES 
AND HCP LATTICES 
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•\/3a 



(53b) 



Another example of a lattice with a basis is the hep 
system, which can be considered_as interpenetrating sim- 
ple hexagonal Bravais latticestj (i.e., stacked triangu- 
lar nets). The primitive vectors of the simple hexag- 
onal Bravais lattice can be taken as ai = ax,a.2 = 



bs = — z 

c 



(53c) 



The reciprocal space is another simple hexagonal lattice, 
with lattice constants in the xy-plane and ^ in the 
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z-direction, rotated by 30° from the real space lattice. 
The Brilluoin zone Wigner-Seitz cell is a hexagonal cylin- 
der, however, for the purpose of the integral needed here 
(Eq. ^) it is more convenient to do the summation in- 
side a cell bounded in the xy-plane by a rhombus formed 
by bi and b2 (See Fig. ^.). The hexagonal cylinder and 
rhombical cylinder cells have equal areas and are equiva- 
lent to each other by appropriate symmetry operations. 
This rhombical cylinder cell is very convenient for evalua- 
tion of the integral /, especially with the variable change 
on q: 



q = xhi + yh2 + zh^ 



(54) 



where the dimensionless parameters x, y, z all range from 
to 1, mapping out the entire cell. This leads to 

(fiq = bi • (b2 X ba) dx dy dz = {2tt)^— dx dy dz (55) 



V 



and the integral is simplified to. 



G^^ terms. Further consideration leads to the nonzero 
coupling elements. 



G^^ = Gf" = J,y,r^ ±ai,±a2,±(ai -aa) (58a) 



^BB 



G'^f = Gf^ = Jz, r = 0,ai,a2,a3,a3 +ai,a3 + a2 

(58b) 

It is notable that it is the first example where the diagonal 
elements are nonzero. The terms where r = are the 
coupling within the basis. 



The q-space couplings (Eq. 37) are found to be 



Gq'^ = Gq ^ = 2Jr,y [cos q • ai + cos q • a2 
-l-cosq • (ai - a2)] 



G^^ = (G^^)* = J4i 



-iq-(a3+ai) , 



-iq-(a3+a2)l 



(59) 



(60) 



dx 



dy 



dz i 1 



{> 



cos 2Trx 



■ cos ^Tiy + cos 27r(x — y) -f A cos 27rz | 



(56) 



The integral gives the correction to the mean-field pre- 
diction, i.e.. 



knTr 



MF 



2(3 + A)J,y52 



(57) 



The integral / was evaluated by the numerical tech- 
niques described above (Sec. for a range of 
anisotropy parameter < A < 2, see Fig. ^. At 
the isotropic limit, A = 1, we get / = 1.44930 and 
nksTc = 5.5199JS'2. In the hmit A ^ 0, the system 
becomes two-dimensional, the integral / diverges loga- 
rithmically due to small-q contributions, and ORF is not 
applicable. 



and what is needed for evaluation of Tci 

G^^G^-^ = Jf [6 -I- 4{cos q • ai + cos q • a2 

-t-cosq • (ai - a2)}] X [l + cosq • as] (61) 



The q=0 exchange strength (Eq. 35) is seen to be 
Jo = G{Jxy + Jz) = 6J2;j,(l -I- A) and determines the 
A-dependent mean-field critical temperature. The spe- 
cific volume per site (for the underlying simple hexagonal 

Bravais lattice) isV/N = ^a?c. The reciprocal space is 
that of the si mple hexagonal Bravais lattice as described 
in Sec. VIA. Therefore, using the rhombical cylinder 



Brilluoin zone cell, the integral / of Eq. ^ can be trans- 
formed using Eqs. p4 and pq into 



1 .1 .1 
dx I dy I dz 

"'0 "'0 



•^0(^0 — Gq^) 



(Jo - G^^)2 - G^BQBA 



(62) 



where we have the transformed quantities 



B. Hexagonal Close Packed Lattices 



G, 



AA 



2Jx 



cos 2'Kx -t- cos 2TTy -\- cos 2'k{x — y) (63a) 



Again there are 6 neighbor displacements from 
an arbitrary site to neighbors in the same plane 
{±ai, ±a2, ±(ai — a2)}, with exchange strength, J^y- 
The difference from the simple hexagonal lattice, is that 
there are 3 neighbors in a layer above and 3 neighbors 
in a layer below the one being considered, with exchange 
couplings Jz: giving 12 neighbors in all. However, to eval- 
uate the matrix elements Gr'-* , we need to consider these 
couplings from the point of view of the simple hexagonal 
Bravais lattice with a basis (see Fig. |[). Thus, an arbi- 
trary A-site, has the 6 neighbor displacements to other 
A-sites in the same xy-plane: {±ai, ±a2, ±(ai — a2)}. 
which will give nonzero G^^ coupling terms. The neigh- 
bors in adjacent planes are B-sites, leading to nonzero 



cos 27rz 

3 -f 2|cos27ra; -I- cos27ry -|- cos27r(a; — y)| (63b) 

The correction to the mean-field prediction is 

6(l + A)J,^g^ 



knTr 



ksT^'^ 



(64) 



/ was evaluated by the numerical techniques described 
above (Sec. VIA), including the 00 extrapolation. 

For example, at the isotropic limit, A = 1, we get / = 
1.34466, and nksTc = 8.92418JS'^ It is interesting to 
note that the same value of / results for the fee lattice 
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(Eq. when evaluated to the same precision. Thus 
the ORF corrections to the mean-field Tc for fee and hep 
lattices, both with 12 nearest neighbors, are the same. 
Some other hep results in a limited range of anisotropy 
A are shown in Fig. ||. Once again, in the limit A ^ 
there is a weak divergence of the integral as the system 
crosses over into a two-dimensional one, with Tc — > over 
a very narrow range of A. ORF is not applicable in this 
limit; Tc should pass ovpito the finite value for the 2D 
triangular lattice model.Ej 

There are a few theoretical results to compare with 
for the hep system. Iij-ajSeries of papers, Domany, Gu- 
bernatis and Auerbachli3~t3 analyzed a Lifshitz tricritical 
point for the hep Ising model, which occurs at a negative 
value of A. As pari-pf their analysis they applied Monte 
Carlo calculationstH to determine the phase diagram; 
very roughly for A = 1 they obtained fc^Tc w 10 J S^. 
Values of Tc at other anisotropics also can be estimated 
from their Fig. 1 but with poor precision. However, it 
does appear that the Onsager results fall below the Monte 
Carlo estimates of Tc, as expected. j— , 

For the hep Heisenberg model, AdleiO estimated Tc 
by a Green's function approach together with a random 
phase approximation. It is surprising to see that at A = 1 
Alder found the correction to the mean-field Tc to be by 
the factor T(-l) = 1.34 ± 0.005, where T(-l) is a cer- 
tain sum over the Brillouin zone. Making a more accurate 
evaluation of T(— 1) by the techniques described here, we 
get T(— 1) = 1.34466, i.e., a value exactly equal to thn 
correction integral / obtained from the ORF procedure.cj 
In fact, it is easy to show that the expression for T(— 1) 
given by Adler is exactly equivalent to our expression 
( |6^ ) for /, including the anisotropic case, A 7^ 1. There- 
fore, the Green's function approach used there is exactly 
equivalent to the ORF procedure presented here; they 
are different approaches to impose the random phase ap- 
proximation. 

Furthermore, in this level of aporoximation, the ques- 
tion posed by Domb and SykesEl and investigated by 
Adler is answered: Tc for fee and hep Ising models are 
the same, even though the hep lattice is more densely 
packed and might be expected to have a higher Tc. Ap- 
parently a more precise procedure is needed to determine 
whether there is a true difference in their critical temper- 
atures. 



VII. REACTION FIELD AND 
THERMODYNAMIC QUANTITIES AT T > Tc 

It is clear that any quantities such as specific heat, 
magnetization, etc, can be evaluated via the RPA Hamil- 
tonian for temperatures away from Tc, provided that the 
reaction field. A, has been determined. Thus we take a 
few sentences to examine how A can be calculated. 

At the critical temperature, the reaction field as deter- 
mined from Eq. (p3) is seen to be 



Ac = A(Tc) = i(Jo - X(7') = ^(^0 - nksTc/S^) (65) 

For higher temperatures, the constraining equation ( |3^ ) 
to determine A is equivalent to 

y /• _fq l"Xo(G^^-2A) _^ 

N Jbz (27r)3 ^1 „ ^^(GAA _ 2A)]' - xgG^^G^^ 

(66) 

Considering the left hand side as a function of A, one can 
apply Newton's method to search for the A at which the 
function passes through 1. This search can be aided by 
the requirement that the denominator of this integrand 
must be positive everywhere in the BZ, including at q = 
0. This leads to the inequality for T > Tc, 

A(r) >l{Jo- Xo') = liJo - nksT/S') < Ac (67) 

Indeed, we have found as well that A(T) < Ac. 

Results for A(T) for the diamond lattice and hep lattice 
(at A = 1) are shown in Fig. |[ The hep lattice, which 
has higher coordination number, also has the stronger 
reaction field at Tc. As the temperature is increased, the 
reaction fields diminish and become of comparable sizes 
for T > 2Tc. It iS|-also expected that the slope relates 
to the specific heat.Q The graph then reasonably demon- 
strates a larger specific heat for hep compared to dia- 
mond near Tc, in contrast to more similar specific heats 
at higher temperatures. 



VIII. CONCLUSIONS 

We have reviewed the standard Onsager reaction field 
approximation for estimating Tc, and have shown how 
it can be extended to apply to a Bravais lattice with a 
basis. The bee lattice was used as a test case because 
it can be calculated either by the standard approach or 
the new method, when considered as sc with a two-atom 
basis. We used the new method to get Tc for diamond 
and anisotropic hep lattice systems, however, it certainly 
can be extended to more complex systems with a greater 
number of atoms per unit cell. For the hep lattice system, 
the ORF procedure used here was found to be exactly 
equivalent to a^Green's function (plus RPA) approach 
used by Adler Jla While it is an approximate method, it 
does give reasonable estimates of Tc and other quantities 
where other methods may be more cumbersome or time- 
consuming to apply. 
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FIG. 1. Evaluation of the integral / of Eq. |5C| for the 
diamond lattice, for different numbers of grid points A^^, 
along each axis. The errors go as l/N^ and extrapolation 
to iV^ ^ oo gives I = 1.79288. 
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FIG. 2. Wigner-Seitz cells (solid line hexagons) for the 
simple hexagonal Bravais lattice reciprocal space, compared 
with the equivalent rhombic cell (dot-dash) used for integrals. 
Segments labeled a,b,c,d are equivalent by symmetry opera- 
tions. 
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FIG. 4. XY projection of some nearest-neighbor bonds 
between A-sites (solid circles) and B-sites (open circles) in an 
hep lattice. Double sohd lines connect A, B sites at the same 
Bravais lattice point. Solid lines show AA bonds (within the 
planes), dotted lines show AB bonds. 
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